DOE/NASA/1 087-2 
NASA CR-191189 


/// - a 7 

04 P 






Three-Dimensional Modeling of Diesel 
Engine Intake Flow, Combustion and 
Emissions-ll 


R.D. Reitz and C .J. Rutland 
University of Wisconsin-Madison 
Madison, Wisconsin 


(NASA-CR-191189) THREE-DIMENSIONAL 
MODELING OF DIESEL ENGINE INTAKE 
September 1993 FLOW, COMBUSTION AND EMISSIONS-2 

Final Report (Wisconsin Univ.) 

64 p 

G3/07 


Prepared for 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
Lewis Research Center 
Under Contract NAG3-1 087 


N94-14448 
Unci as 
0187798 


for 

U.S. DEPARTMENT OF ENERGY 
Office of Transportation Technologies 





Three-Dimensional Modeling of Diesel Engine 
Intake Flow, Combustion and Emissions - II 

DOE/NASA-Lewis 
University of Wisconsin Grant* 

December, 1992 

Principal Investigators: R.D. Reitz and C.J. Rutland 
Engine Research Center 
Department of Mechanical Engineering 
University of Wisconsin-Madison 

’Research sponsored by the U.S. Department of Energy, Assistant Secretary of 
Conservation and Renewable Energy, Office of Transportation Technologies, administered by 
NASA-Lewis under Grant number NAG 3-1087, William Wintucky Grant Monitor. 

Report Period 10/91 - 12/92 


1 


ABSTRACT 


A three-dimensional computer code, KIVA, is being modified to include 
state-of-the-art submodels for diesel engine flow and combustion. Improved 
and/or new submodels which have already been implemented and 
previously reported are: wall heat transfer with unsteadiness and 

compressibility, laminar-turbulent characteristic time combustion with 
unburned HC and Zeldo’vich NOx, and spray/ wall impingement with 
rebounding and sliding drops. Progress on the implementation of improved 
spray drop drag and drop breakup models, the formulation and testing of a 
multistep kinetics ignition model and preliminary soot modeling results are 
described in this report. In addition, the use of a block structured version of 
KIVA to model the intake flow process is described. A grid generation 
scheme has been developed for modeling realistic (complex) engine 
geometries, and computations have been made of intake flow in the ports 
and combustion chamber of a two-intake-valve engine. The research also 
involves the use of the code to assess the effects of subprocesses on diesel 
engine performance. The accuracy of the predictions is being tested by 
comparisons with engine experiments. To date, comparisons have been 
made with measured engine cylinder pressure, temperature and heat flux 
data, and the model results are in good agreement with the experiments. 
Work is in progress that will allow validation of in-cylinder flow and soot 
formation predictions. An engine test facility is described that is being used to 
provide the needed validation data. Test results have been obtained showing 
the effect of injection rate and split injections on engine performance and 
emissions. 
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INTRODUCTION 


The diesel engine is the leading heavy-duty power plant because of its 
superior energy efficiency. However, because of environmental concerns, 
proposed federal emission standards require reductions in both nitric oxides 
(NOx) and particulates. A detailed understanding of combustion is required 
in order to work effectively at reducing these by-products of combustion 
within the engine cylinder, while still not compromising engine fuel 
economy. 

The objective of this research program is to develop an analytic design tool 
for use by the industry to predict engine performance and emissions. It is 
expected that the use of advanced modeling tools will enable engine 
development times and costs to be reduced. The three-dimensional 
computer code, KIVA [1,2] is being used since it is the most developed of 
available codes. Part of the research involves implementing state-of-the-art 
submodels into KIVA for spray atomization, drop breakup/ coalescence, 
multi-component fuel vaporization, spray /wall interaction, ignition and 
combustion, wall heat transfer, unburned HC and NOx formation, soot and 
radiation, and modeling the intake flow process. 

Previous progress has been described by Reitz et al. [3,4], and in a previous 
Annual Report [5], where it was shown: that adding the effects of 
unsteadiness and compressibility to the heat transfer submodel improves the 
accuracy of heat transfer predictions; that drop rebound should be included in 
spray /wall interaction models since rebound can occur from walls 
particularly at low impingement velocities (e.g., in cold-starting); that the 
influence of vaporization on the atomization process should be considered 
since larger spray drops are formed at the nozzle with vaporization; that a 
laminar-and-turbulent characteristic time combustion model has the 
flexibility to match measured engine combustion data over a wide range of 
operating conditions; and, finally, that the characteristic time combustion 
model can also be extended to allow predictions of ignition. 

The implementation of improved models for the prediction of spray drop 
drag and breakup, and for modeling ignition and combustion, intake flow 
and soot formation is described in this report. The research also involves the 
use of the computer code to assess the effects of subprocesses on diesel engine 
performance, and the accuracy of the predictions is being tested by 
comparisons with engine experiments. To date, comparisons have been 
made with measured engine cylinder pressure, temperature and heat flux 
data, and the model results are in good agreement with experiments [3, 4, 5]. 

Work is also in progress that will allow validation of in-cylinder flow and 
soot formation predictions. An engine test facility is described that has been 
constructed and is being used to provide the validation data. 
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Close contact is being maintained with engine industries during the course of 
the research. This includes participation in DOE diesel engine working group 
meetings 1 . In addition, technology has been transferred directly to industry 
(e.g., we have already made FORTRAN code for the ignition subroutine listed 
in Appendix 3 of this report available to the Software Development Group, 
Engine Division, Caterpillar Inc. - W. Brown of this group spent the past year 
at Madison). Two of our graduate students (D. Nehmer and R. Hessel) are 
now employed by Cummins Engine company and have taken elements of the 
technology with them to Cummins. We have continued our involvement 
with the KTVA Users Group (including co-organizing the 2nd International 
KIVA Users Group Meeting in Detroit on February 23rd, 1992, and producing 
Newsletters 5, 6 and 7, see Appendix 1). This activity has kept us informed of 
code enhancements and new developments elsewhere in the international 
user community. 

BACKGROUND OF RESEARCH PROGRAM 

The research program was initiated in the Fall of 1989 as a five year program. 
The emphasis of the research is on the application of a comprehensive engine 
combustion code to assess the effect of the interacting subprocesses on diesel 
engine performance, rather than on the development of new models for the 
subprocesses themselves. The elements of the code are being assembled from 
existing state-of-the-art submodels. This use of multidimensional modeling 
as an engine development tool is timely and justifiable due to recent 
advances in submodel formulations. 

The computer model is based on the KIVA-2 and KTVA-3 computer codes 
that were written at the Los Alamos National Laboratories [1,2] and allow 
computations of turbulent, multi-phase, combusting 3-D engine flows. The 
codes solve the conservation equations for the transient dynamics of 
vaporizing fuel sprays interacting with flowing multicomponent gases which 
are undergoing mixing, ignition, chemical reactions, and heat transfer. The 
codes have the ability to calculate flows in engine cylinders with arbitrary 
shaped piston geometries, and include the effects of turbulence and wall heat 
transfer. 

Many of the in-cylinder processes occur on time and length scales that are too 
short to be resolved in practical computations. Thus, submodels are needed 


!prof. Rutland gave the presentation "Diesel Intake Flow Modeling," Spring 
1992 DOE Meeting, Southwest Research Institute, San Antonio, Texas, May 4-5, 
1992. Prof. Reitz gave the presentations "Progress in 3-D Modeling of Diesel 
Engine Intake Flow, Combustion and Emissions' and "Measurements of the 
Effect of Injection Rate and Split Injections on Diesel Engine Soot and NOx 
Emissions" at the Fall 1992 DOE Diesel Group Meeting, Princeton University, 
November 5-6, 1992. 
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to describe the processes, and the submodels require detailed validations 
before they can be used for engine performance predictions (Reitz [6]). For 
example, for modeling sprays, submodels are needed to describe drop drag, 
drop turbulence interaction, vaporization, atomization, breakup and 
coalescence, and spray/wall interaction. These spray processes influence the 
combustion and pollutant formation processes profoundly through their 
effect on the fuel-air distribution in the engine. 

Individual submodels are usually developed starting from simplified 
theoretical models or experiments that isolate a particular process. However, 
the accuracy of predictions made using combinations of submodels to describe 
the performance of the overall engine must be assessed by comparisons with 
detailed and informative engine experiments. An important component of 
the present research is to use engine comparisons to identify areas where 
submodels need improvement, and to develop improved models where 
necessary. 

The program was expanded to include the modeling of the diesel engine 
intake process in 1990 due to its importance in the combustion process. In 
1991 the research program was expanded again to include engine experiments 
for the acquisition of validation data. A summary of the 1989-1991 research 
activities has been described in the 1991 Annual Report, Ref. [5]. 

The research program has been organized into three main activities: Part A - 
Submodel Implementation, Part B - Model Application and Part C - 
Experiments. Progress in the Part A tasks is described on page 7 and in Table 1 
below. Part B tasks are described on page 27. The Part C tasks are described on 
page 29. 


Table 1 Status of K3VA Submodels 



Original KIVA 

Submodels being 
implemented / upd ated 

F 

combustion 

Arrhenius 

laminar-turbulent char time 

1-3 

unburned HC 

Arrhenius 

laminar-turbulent char time 

1 

NOx 

Zeldo'vich 

Zeldo'vich 

1 

heat transfer 

law-of-the-wall 

compressible, unsteady 

1-2 

wall impinge 

none 

rebound-slide model 

1-2 

atomization 

Taylor Analogy Breakup 

surface-wave-growth 

1-3 

intake flow 

assumed initial flow 

KIVA-3 computed intake 
flow 

2-4 

ignition 

none 

Shell multistep model 

2-3 

vaporization 

single component 

multicomponent 

1-4 

soot 

none 

KIVA-coal 

3-4 

radiation 

none 

Rayleigh-limit, band model 

3-4 
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CURRENT PROGRESS 


PART A - Submodel Implementation 

A summary of the submodels that are being added to KIVA and the time- 
table for their implementation is given in Table 1. The models include 
improved models for combustion, unburned HC and NOx, wall heat transfer, 
spray /wall impingement, atomization, intake flow, ignition, multi- 
component vaporization, and soot and radiation. Details of the current and 
projected status of the models are presented in the last column of Table 1. 

The research program is currently in its fourth year. Progress to date in each 
of the Part A tasks scheduled for the third year of the grant is discussed in this 
section. Additional details of progress are given in Ref. [5]. Appendix 2 gives 
additional details of published papers connected with the research activities 
in the past year. Appendix 3 contains a FORTRAN source code listing for the 
multistep ignition submodel that has been implemented as part of the 
research. 

Improved and/or new submodels which have already been implemented in 
KIVA have been described by Reitz et al. [3,4], including models for: wall heat 
transfer with unsteadiness and compressibility, laminar-turbulent 
characteristic time combustion with unburned HC and Zeldo'vich NOx, and 
spray/wall impingement with rebounding and sliding drops. In addition, 
progress on intake flow modeling, including the generation of grids for 
computations in realistic manifold geometries, and the use of an block 
structured mesh version of KTVA, KTVA-3, have been described by Reitz et al. 
[4]. 

Results to date show that adding the effects of unsteadiness and 
compressibility improves the accuracy of the heat transfer predictions; spray 
drop rebound can occur from walls at low impingement velocities (e.g., in 
cold-starting); and that a laminar-and-turbulent characteristic time 
combustion model has the flexibility to match engine combustion data over a 
wide range of operating conditions. Continued progress on intake flow 
modeling, and on diesel spray atomization and ignition/combustion and 
emissions models is described next. 
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ATOMTZATION AND DROP DRAG - As noted by Reitz et al. [4], vaporizing 
spray computations have been found to be very sensitive to numerical 
parameters, especially to the grid size near the exit of nozzle. This motivated 
re-examination of the parameters of the atomization model since the 
atomization process occurs mostly near the nozzle exit. For this purpose, 
comparisons have been made by Liu et al. [7] with drop breakup, velocity and 
trajectory measurements of Liu and Reitz [8]. 

The experiments consisted of a Berglund-Liu drop generator and an air jet 
nozzle, arranged in a cross-flow pattern so that the drops are suddenly 
exposed to the air as in Fig. 1. Drop sizes were measured with a Phase 
Doppler Particle Analyzer [8]. This experimental configuration is useful 
because the relatively low liquid and spray drop density does not favor 
collisions and coalescences of the drops which complicate the interpretation 
of typical diesel spray experiments. Drop breakup was. modeled using both 
the wave instability model (Reitz [9]) that has benn used in our diesel engine 
combustion simulations, and the standard KIVA TAB breakup model [10]. 
The breakup models contain model constants that were also evaluated as part 
of the present study. 

The computations shown in Fig. 1 were made on a 3-D mesh of 32x16x84 cells 
in the radial, azimuthal and axial directions, respectively. Figure 1 shows 
computed drop locations and gas velocity vectors in the plane of the nozzle, 4 
ms after the start of the injection for a case with an air velocity of 100 m/s at 
the air nozzle exit. 



Fig. 1 Computed drop locations and gas velocity vectors 4 ms after the start of 
injection for air jet velocity of 100 m/s. Stream of 170 |im diameter drops 
enters air jet from the left at 16 m/s. 
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Measured and predicted (parent or outermost) drop trajectories are compared 
in Fig. 2. In this figure the drops enter from the left (at X=Y=0), and the 
measured data is indicated by squares. Figure 2 shows that good agreement 
with the measured trajectories can be obtained when a modified drag 
coefficient is used that accounts for the fact that the parent drops are deformed 
by the gas flow as they breakup. The drop deformation is modeled 
dynamically as a function of the local flow conditions around each drop [7]. 

The corresponding predicted drop size variation across the air jet using the 
two different breakup models is shown in Fig. 3. In the wave instability 
model (Wave [9]), drop breakup times and breakup sizes are assumed to be 
proportional to wave growth rates and wavelengths obtained from a stability 
analysis of liquid jets. The final (stable) drop size reached by drops that 
penetrate all the way through and emerge from the other side of the air jet are 
seen to be in good agreement with the measured drop size at the jet edge in 
Fig. 3. Figure 3 also shows computations made with the standard KIVA TAB 
breakup model [10]. The TAB model final drop size prediction is seen to 
under estimate the measured drop size significantly. 


case 2 



Fig. 2 Comparison of wave breakup model predictions and measured drop 
trajectory for air jet velocity 59 m/s. Solid line - standard drop drag, dashed 
line - dynamic drag model that accounts for drop distortion effects. 
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case 2 



Fig. 3 Predicted drop Sauter mean diameter variation with distance across the 
jet for air jet velocity 59 m/s using the wave (solid line) and TAB (dashed 
line) models. Solid circle shows PDPA measured drop diameter outside air 
jet. 

The effect of the sensitivity of drop trajectories to drop drag coefficients 
observed in Fig. 2 for single drops was also assessed for diesel-type sprays [7]. 
The results show that spray-tip penetration is quite insensitive to the value of 
the drop drag coefficient. This somewhat surprising result was found to be 
due to the fact that changes in the drag coefficient produce changes in the 
drop-gas relative velocity which, in turn, cause changes in the spray drop size 
through the drop breakup and coalescence processes. The changes occur in 
such a way that the net effect on the spray tip penetration is small over the 
range of tested conditions. However, the results indicate that the drag model 
does influence the spray drop sizes, and this has an important effect on fuel- 
air mixing since spray vaporization rates depend strongly on the drop size. 
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TCNITION AND COMBUSTION MODELS - The ignition delay is an 
important parameter in the operation of diesel engines since it influences 
hydrocarbon and NO* emissions. During the delay period, the injected fuel 
undergoes complex physical and chemical processes such as atomization, 
evaporation, mixing and preliminary chemical reaction. Ignition takes place 
after the preparation, and reaction of the fuel-air mixture leads to fast 
exothermic reaction. Approaches to describe autoignition phenomena in 
diesel engines in multi-dimensional models have been reviewed by Kong 
and Reitz [11]. 

The kinetics chemistry submodel used for ignition and combustion modeling 
in the standard KIVA considers a single step Arrhenius mechanism for the 
stoichiometric reaction of the fuel. Kong et al. [12] used the Arrhenius model 
to describe ignition together with a laminar and turbulent characteristic time 
model for combustion. The latter model has been used successfully in spark- 
ignited engine studies and for premixed- and fuel-injected two-stroke engines 
(Reitz [13] and Kuo and Reitz [14, 15]). In this way the combustion model has 
been extended to allow predictions of ignition. Kong et al. [12] found that it 
was possible to match homogeneous-charge, compression-ignition engine 
experiments reasonably well with one set of model constants. 

A more accurate prediction of ignition has been achieved by implementing 
the multistep Shell kinetics model of Halstead et al. [16, 17] in KTVA. The 
implementation details of the model are given in Appendix 3 of this report. 
Figure 4 shows a comparison between spray ignition experiments and KIVA 
computations of ignition using the Shell model. The experiments were 
conducted in a constant volume bomb at the Sandia Labs (Edwards et al. [18]), 
and used diesel fuel (simulated in the computations as Hexadecane) injected 
into compressed air at 30 atmospheres and 1070 K. Ignition is seen in Fig. 4 to 
have occurred by 1.13 ms after the beginning of the injection, as represented 
by a luminous zone located at the edge of the spray, some distance 
downstream of the nozzle. By 1.51 ms the luminous zone extends over most 
of the spray, indicating that combustion has commenced. 

The predicted temperature contours are also shown in Fig. 4 for comparison, 
and it is seen that the location and timing of the ignition process is predicted 
fairly well by the model. Similar good levels of agreement between measured 
and predicted ignition delay times were also found for other chamber gas 
temperatures as shown in Fig. 5. 

Figure 6 shows results of comparisons of measured engine cylinder pressure 
data of Yan and Borman [19] in a Cummins NH engine operated under the 
conditions given in Table 2. The predictions were made using the standard 
Arrhenius and the multistep Shell ignition models. Post-ignition 
combustion was simulated using the single-step Arrhenius model in both 
cases. The single-step kinetics model is seen to work reasonably well. 



L13ms 



1.51 ms 


L=641 K, H=2320 K 



L=664K, H=2390K 


Fig 4. Comparison of measured and predicted spray ignition data for the 
Sandia constant volume bomb experiments [17]. 


However, it was found necessary to use a different combustion (pre- 
exponential) model constant than that used for computations with the 
Caterpillar engine (see later Fig. 15), indicating that the single-step Arrhenius 
model is not predictive. 
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1000/T 
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Core Temperature at injection (K) 


Fig 5. Comparison of measured and predicted ignition delay times for 
ignition for the Sandia constant volume bomb at different chamber gas 
temperatures. 


Table 2 Cummins Engine 
Cylinder Bore 
Stroke 

Compression ratio 
Displacement 

Number of spray nozzle orifices 

Nozzle hole diameter 

Engine speed 

Overall equivalence ratio 

Air flow rate 

Fuel flow rate 

Intake tank pressure 

Intake tank temperature 

Cylinder wall temp 

Cylinder head 

Piston surface 

Mass average gas temperature at IVC 
Cylinder pressure at IVC 
Swirl number 
Fuel 

Injection starts 
Injection ends 


Data 

139.7 mm 
152.4 mm 
13.23 

2.33 liters 
8 

0.2 mm 

1500 r/min (constant) 
0.6 

0.353e-2 kg/cyde 
0.144e-3 kg/cyde 
148.2 kPa 
302 K 
405 K 
486 K 
578 K 
359 K 
157.9 kPa 
1.0 

Tetradecane 
18° BTDC 
11° ATDC 
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Fig. 6 Measured and predicted pressure histories with the standard one-step 
Arrhenius model (Case 1) and the multistep Shell ignition model (Case 2). 
Cummins engine at 1500 rev/min. 

The multistep ignition model is more predictive. The computations of Fig. 6 
used parameters derived from out-of-engine, low temperature chemistry 
data, and the multistep model is seen to give excellent results prior to the 
main stage of combustion. 

Details of the predicted location and time-history of the ignition process are 
summarized in Fig. 7. Fuel injection starts at 18 degrees BTDC. Ignition takes 
place some distance downstream of the nozzle at the spray edge at about 9 
degrees BTDC. This is evidenced by the high temperature contours which are 
seen to surround the spray in this region in Fig. 7, consistent with the 
constant volume bomb experimental ignition results of Fig. 4. The build-up 
and subsequent rapid consumption of the intermediate, branching and radical 
species that control the ignition process are shown as a function of crank 
angle in the ignition cell in the inset plot. After the ignition, the gas 
temperature variation in the ignition cell reflects the balance between the 
energy released due to combustion and the energy required to vaporize the 
liquid fuel. 
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Fig. 7 Predicted spray drop locations and temperature contours in the 
Cummins engine at 8 degrees BTDC. Inset plot shows variation of gas 
temperature and ignition species concentrations with crank angle in the 
ignition cell. 

The discrepancy between the predictions of both the standard Arrhenius and 
the multistep kinetics models and the measured results after ignition suggests 
the need to consider other effects such as turbulent mixing on post-ignition 
combustion. Also, discrepancies with measured data have been found to be 
influenced by numerical errors, and by the details of the fuel injection 
parameters by Reitz et al. [4]. Further investigation of these effects is in 
progress. 
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SOOT AND NOX MODELS - As described in the previous annual report, the 
extended Zeldo'vich NOx model [20] has also been implemented in KIVA. 
Two soot formation models have been implemented into KIVA in the 
current activity. The first is the Surovkin model [21] which had previously 
been implemented by Gentry et al. [22] in the KIVA-COAL code (for coal- 
slurry combustion modeling). The second is the simpler Hiroyasu model 
which is a two equation model for the soot mass formation rate and the soot 
mass oxidation rate and has been implemented in KTVA by Bertoli et al. [23]. 

In the Surovkin model the soot formation process starts by considering the 
formation of radical nuclei which grow to become soot particles by surface 
growth and agglomeration mechanisms. The model consists of 5 
simultaneous rate equations for radicals, radical size, particles, particle size, 
and a hydrocarbon consumption term. The soot particles themselves are 
subject to oxidation processes following the Nagle and Strickland-Constable 
mechanism [24]. In the present implementation the equations were solved 
using the Runge-Kutta approach, since the original explicit method of Ref. 
[22] was found not to work in the high fuel concentration and high 
temperature environment of the diesel engine. 

The results given in Fig. 8 show predictions of the thermal decomposition of 
C6H6 in an inert high temperature environment, and the results agree well 
with the model results of Surovkin [21] (not shown). However, detailed 
investigations of the results suggest that the soot particle size model needs 
further development, especially for engine applications, since the model can 
produce more soot mass than the original fuel has available. Although the 
Surovkin model has been used in engine computations (see Fig. 9), it is 
believed that the model over-predicts the soot mass. One possible solution 
would be to modify the model constants. 

Preliminary calculations with the Hiroyasu model indicate that this model 
computes soot mass quantities which are in the known engine range. 
Although this model was not developed for modeling the thermal 
decomposition of a hydrocarbon fuel, results obtained using the model to 
compute the soot mass for the thermal decomposition of C6H6 in the absence 
of oxygen are compared with the Surovkin model in Fig. 8. The Hiroyasu 
model predicts soot mass levels that are about two orders of magnitude lower 
than those from the Surovkin Model. 

Thus it is expected that the Hiroyasu model will work well in the engine 
environment, but the Surovkin model, originally developed for thermal 
decomposition, will require additional work to be useful for engine 
computations. 

Nevertheless, the results in Fig. 9 show predicted temperature, NOx and soot 
particle distribution contours obtained using the Surovkin model in the 
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Cummins engine at 6 degrees BTDC. The temperature data indicates that a 
combustion region surrounds the spray near the nozzle ('h' identifies the 
high contour location). The NOx contours reveal that the NOx 
concentrations are highest in region where the gas temperatures are the 
highest, as expected from the sensitivity of the NOx reaction mechanism to 
temperature. The soot particle size contours indicate that the soot 
concentrations are largest in the colder regions near the spray tip. It is worth 
noting that soot and NOx are predicted to be formed in different regions of 
the combustion chamber. This observation is important in developing an 
understanding of the well-known soot-NOx trade-off. 



Fig. 8 Comparison of Hiroyasu and Surovkin soot models for the high 
temperature thermal decomposition of C6H6 in the absence of oxygen. 
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INTAKE FLOW MODELING - Significant advances have been made in intake 
flow modeling. A methodology has been established by Rutland et al. [25] for 
generating a realistic grid for the complex geometry of the intake port on the 
Caterpillar head as shown in Fig. 10. 

Creation of the grid required a digitized representation of the geometry. Next, 
the geometry was imported into the GRIDGEN code of Steinbrenner et al. [26]. 
Several tools were developed as part of the grid creation process, including a 
reshaping technique where two or more initial grids are created, each with 
some aspects desired in the final grid. The result is a composite grid that has 
the required characteristics, and is much simpler than what could be achieved 
with GRIDGEN alone [25]. A grid management code was created to help 
transfer the GRIDGEN grids to KIVA-3. The reshaping techniques are 
implemented in this program along with the ability to manipulate grid blocks 
and boundary conditions. The output from this program is run through a 
modified version of the KIVA grid generation program, K3PREP. In K3PREP 
the grid generation is disabled but the cell and vertex flags required by KIVA 
are initialized. 

Computations have been made for steady flow (fixed valve lift) in the port, 
across the two-intake valves and into the cylinder. Pathlines showing 
predictions of the mixing and flow structures that are setup in the cylinder are 
shown in Fig. 11. The pathlines are obtained using 'particle' tracks in an 
instantaneous velocity field. It can be seen that the flow structures change 
through the cylinder and are not readily described by single parameters such 
as swirl or tumble ratios. 

The steady intake flow computational results have also been compared with 
experimental mass flow measurements obtained using a SuperFlow 600 flow 
bench. The calculations were made by specifying the inlet and outlet 
pressures and calculating the mass flow rates at steady state. The results in 
Table 3 show that there is reasonably good agreement between the calculated 
and experimental mass flow rate results. 

Table 3. Comparison of mass flow rates from experimental flow bench 

and KIVA-3 results. 


Case 

Experiment* 

KIVA-3 

% Difference 

1 

37.68** 

40.69 

-7.99 

2 

54.17 

52.00 

+4.01 

3 

45.96 

45.24 

+1.57 

4 

62.33 

57.16 

+8.29 


* mass flow rate in(g/sec) 

**accuracy estimated as +/- 8% due to leakage 
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Fig. 10 Computational grid used for the two-intake valve Caterpillar engine. 
The bell-mouthed entrance to the intake manifold is at the intake surge tank. 
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Pig. ii Predicted in-cylinder fluid pathlines for steady intake flow in the 
Caterpillar engine. 
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The capability to accurately simulate moving valves has also recently been 
implemented. This work takes advantage of the block structured gridding 
and the 'snapper' algorithm of KTVA-3 [2]. The implementation of the 
snapper algorithm was described previously by Reitz et al. [4] and it permits 
regions of the grid to move relative to each other. Grid point connectivity is 
maintained by periodically regridding or snapping. Adapting the algorithm 
for use around the valves has proven to be very successful. It is 
computationally efficient, accurate, and there is no distortion of the physical 
surfaces. 

In addition, a method for closing the valves has been developed. The contact 
between the valve and the seat is a region rather than a point. This can result 
in severe time step limitations and difficulty in dealing with the trapped mass 
since the cell volumes go to zero. The problem has been overcome by using 
walls to block the flow. The walls are erected when the flow rate across the 
valves is nearly zero, but before the valve actually touches the valve seat. 
The inaccuracies introduced are negligible because the flow rates are so low. 
This approach avoids problems with time step limitations and mass 
conservation. 

Computations have been made with dual moving and closing valves. The 
valves move according to individual, specified lift profiles. A constant 
pressure boundary condition was used at the top of the port at the intake 
surge tank. Results from a sample calculation are shown in Fig. 12 and the 
run conditions are shown in Table 4. Fig. 12 shows the velocity field on a 
vertical diametral plane that passes through one of the valves. 

A piston with a bowl was added to the grid as shown in Fig. 12. This required 
some modification of the combustion chamber grid so that it would conform 
to the shape of the bowl. This modification also required some changes to the 
grid around the valve curtain region. These changes were implemented 
fairly rapidly using the capability of GRIDGEN to incorporate an existing grid 
as a geometry database. 

The intake valves open at 8.5 degrees before top dead center. The velocity 
vector results in Fig. 12 at 0 degrees (TDC) show that fluid flows back into the 
intake manifold initially as the valves begin to open. By 60 degrees ATDC the 
piston is moving downward and strong intake flow jets across the valves 
create large scale vortical flow patterns in the engine. 
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Fig. 12 Velocity vectors in a diametrical plane through one intake valve 
during the intake stroke for dual moving valves implemented in KTVA-3; (a) 
0 degrees (TDC) (b) 60 degrees ATDC Caterpillar engine 1600 rev/min. 


Table 4 Parameters for the unsteady intake flow calculations 


Engine R.P.M. 

1600 

Calculation crank angles 

10° BTDC to 3350 ATDC 

Intake valve opening 

8.5° BTDC 

Intake valve closing 

184° ATDC 

Initial pressure (port and chamber) 

179 KPa 

Initial temperature 

723 K 

Initial turbulent kinetic energy 

3.7m 2 /s 2 

Boost pressure at inlet boundary 

184 KPa 
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Volume averaged flow quantities as a function of crank angle are presented 
in Fig. 13. The turbulent kinetic energy is seen in Fig. 13a to increase as the 
intake flow jets create regions of high production. However, the total 
turbulence levels begin to decrease even before the intake valves close. By 
the time of fuel injection the average turbulence level has decayed by an 
order of magnitude from the peak value. 

The volume averaged turbulence integral length scale also increase during 
the intake process as shown in Fig. 13b. However, this lags the turbulent 
kinetic energy and peaks around the time of intake valve closing. This lag is 
consistent with turbulence being created at smaller scales (in the shear layers 
of the intake jets) and decaying into larger scales as production diminishes. 
The turbulence length scale also decreases significantly by the time of 
injection. 

The volume averaged swirl ratio is presented in Fig. 13c. The swirl ratio is the 
average angular velocity of the fluid normalized by the engine r.p.m. with 
negative values indicating clockwise rotation. The magnitude of swirl 
increases during intake, decreases slightly before intake valve closing and 
remains nearly constant for the rest of the compression calculation. 

Since the steady flow calculations indicated that a global swirl parameter does 
not provide an adequate means of assessing intake flow structures, more 
detailed information about the swirl is of interest. Figure 14 shows swirl 
profiles averaged over horizontal planes, and plotted as a function of vertical 
position in the chamber and crank angle. In these figures the piston face is 
indicated by an arrow. Non-zero swirl below the piston face occurs in the 
piston bowl. 

The plots in Fig. 14 show the rapid increase of swirl during intake. However, 
most of the increase is located close to the piston face. This is consistent with 
the findings of the steady flow calculations (Fig. 11). The flow from the intake 
jets creates localized vortical structures and prevents a chamber-wide 
structure from forming near the head. Close to the piston face, there is no 
interference from the valve flow and the flow is constrained to lie in a 
horizontal plane. 

The vertical distribution of swirl can be compared to the volume averaged 
values in Fig. 13c to show how it changes with time. This is indicated in 
Table 5 in which the peak swirl magnitudes are compared to the volume 
averaged magnitudes. The ratio of peak to mean value diminishes from 3.6 
near the end of intake valve closing to 1.5 at start of injection. The large ratio 
near intake value closing is due in part to swirl in the upper part of the 
chamber that is in an opposite direction from the predominate swirl direction 
lower in the chamber. The magnitude of this opposite swirl is not large but it 
occurs over a large region and results in a lower mean magnitude. As the 
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Crank Angle, degrees 


(C) 

Fig. 13 Volume averaged turbulence intensity (a), length scale (b) and swirl 
ratio (c) as a function of crank angle during the intake stroke. 



(a) (b) 

Fig. 14 Swirl ratio variation with depth below the head at (a) 60 and (b) 180 
degrees ATDC. The arrow indicates the piston face. Non-zero swirl ratios 
below the face are in the piston bowl. 
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Table 5 Swirl magnitudes during intake. 


Crank Angle 

Peak Swirl 
Magnitude 

Mean Swirl 
Magnitude 

Ratio 

180 

2.5 

0.7 

3.6 


1.4 

0.6 

2.3 

335 

0.8 

0.55 

1.5 


cycle continues this opposite swirl diminishes until it is completely gone by 
the time of injection. 

The above results indicate that engines flow fields contain fairly complex 
structures. In general, the flows were regular and smooth upstream of the 
valves. Downstream of the valves in the cylinder, the flows progressed from 
complicated regions with smaller scales to more regular regions with 
prominent swirl patterns. Simple, global parameters such as swirl or tumble 
do not adequately describe these flows. The vertical structure of the planar 
averaged swirl shown in Fig. 14 provides insight into evolution of swirl 
during the intake and compression processes. 
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PART B - Model Application 

PERFORMANCE TESTS - A preliminary set of performance tests have been 
conducted on the Caterpillar engine (see PART C - Engine Experiments), and 
the results have been compared with KIVA predictions. The experiments 
used the standard Caterpillar pump-line injection system and the 
experimental conditions were chosen to match data obtained from Caterpillar 
on a multi-cylinder version of the engine. The experimental conditions are 
given in Table 6. The results were found to be in good agreement with the 
multi-cylinder engine data when allowance was made for the increased 
friction work of the single-cylinder engine. 

A diesel engine simulation program was also used to estimate the cylinder 
wall, valve and piston temperatures, and the simulation was modified to 
account for the engine geometry, intake and exhaust surge tank volumes, and 
its heat release model constants were optimized to match the measured in- 
cylinder pressure trace. 


Table 6 Caterpillar Engine Data 


Cylinder Bore 137.19 mm 

Stroke 165.1 mm 

Compression ratio 15.0 

Displacement 2.44 liters 

Number of spray nozzle orifices 6 


0.1295 mm 

1000 rev /min (constant) 
0.564 

0.267e-2 kg/cycle 
0.883e-3 kg/cyde 
124.0 kPa 
323 K 
422 K 
586 K 
578 K 
773 K 

Mass average gas temperature at IVC 366 K 


Cylinder pressure at IVC 133.0 kPa 

Swirl number 0.978 

Fuel Tetradecane 

Injection starts 23 s BTDC 

Injection ends 6.85° ATDC 


Nozzle hole diameter 
Engine speed 
Overall equivalence ratio 
Air flow rate 
Fuel flow rate 
Intake tank pressure 
Intake tank temperature 
Cylinder wall temp 
Cylinder head 
Piston surface 
Valves face 
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The corresponding predicted and measured cylinder pressure versus crank 
angle data are shown in Fig. 15. As can be seen, there is good agreement 
between the experiments and the simulation, and the KIVA computed 
pressures. It should be noted Aat the KIVA results were found to be sensitive 
to grid resolution and the injection details, and the results shown in Fig. 15 
used the fine grids recommended by Gonzalez et al. [27], and measured 
injection duration data. 

Work is currently in progress on assessing the accuracy of KIVA predictions 
of the effect of injection rate and split injections using the experimental data 
to be described in the next section. 


Predicted and Measured Cylinder Pressure 
Caterpillar Engine 1000 rev/min 



Fig. 15 Comparison between KIVA and (zero-dimensional) engine 
simulation predictions and the experimental cylinder pressure versus time 
for the Caterpillar engine at 1000 rev/min. 
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PART C - Experiments 

A supporting experimental program has been initiated to generate data 
needed for KIVA validations. The engine experiments are conducted using a 
single cylinder Caterpillar 3406 oil- test research engine. In addition to 
standard measurements of cylinder gas and fuel pressures, and ot er 
measurements, combustion visualization and photography experiments are 
being made through a window that replaces one of the exhaust valves. In 
setting up the experiments, effort has been made to ensure that the 
modifications to the engine for optical access are minimal so that the results 
will be representative of the actual engine. 

The engine has also been equipped with a Nippondenso U2, electronically 
controlled ultra-high-pressure fuel injector with programmable injection rate 
shapes and timings, and injection rate data is obtained with a Bosch injection- 

rate meter. 

The experiments are being conducted on a timetable that ensures that data 
will be available coincident with the model implementation phases of Table 
1. Details of progress on the experiments are summarized below. 

ENGINE FACILITY - A photograph showing the Caterpillar single cylinder 
research engine is presented in Fig. 16. An electrical air heater supplies 
heated, pressurized air (to simulate turbo-charging) to a large volume intake 
surge tank (partially visible at upper right in Fig. 16). A short-length intake 
duct with a bell-mouthed entrance leads from the bottom of the intake duct to 
the engine cylinder head. The exhaust port leads to the insulated exhaust 
settling chamber (upper left in Fig. 16). Details of the engine geometry are 
summarized in Table 7. 
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Table 7 Caterpillar Engine Specifications 


Bore 

Stroke 

Connecting Rod Length 
Displacement 
Compression Ratio 
Piston crown 
Engine Speed 
Intake pressure 
Intake temperature 
Intake valve close timing 
Swirl ratio (nominal) 


137.19 mm 
165.1 mm 
261.62 mm 
2.44 L 
15 

Mexican hat 
up to 2100rev/min 
up to 250 kPa 

up to 450 K 
147 deg BTDC 
0.978 



Fig. 17 Optical arrangement for optoacoustic gas temperature measurements. 

In Fig. 16, the engine is shown with the cyl inder hea d water cooling removed 
for motoring engine optical exper iments. The Nd:Yag laser (visible at the 
front left) is being used for in-cylinder velocity and temperature 
measurements with optical access through a quartz window that takes the 
place of one of the exhaust valves as indicated in Fig. 17. 

A new technique which is based on optoacoustic phenomena has been 
developed for measuring in-cylinder gas temperature and turbulent 
diffusivity and is described by Giangregorio et al. [28]. The pulsed laser beam 
is focused to cause local ionization of air at a point in the combustion 
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chamber as shown in Fig. 17. This initiates a shock wave and creates a hot 
spot. The local temperature and turbulent diffusivity are determined by 
monitoring the shock propagation and the hot spot growth, respectively, with 
a schlieren photography system. 

Measured gas temperatures have been compared with KIVA predictions as 
shown in Fig. 18 for a motored engine case at 750 rev/min. As can be seen, 
the agreement is good. The extension of the method for the measurement of 
turbulence diffusivity is currently in progress. 

EMISSIONS MEASUREMENTS - A full dilution emission tunnel has been 
constructed for the engine experiments. The tunnel allows soot emissions 
(dry and SOF using Soxhlet extraction) from the engine to be studied as a 
function of the engine operating conditions. Gaseous emissions are also 
monitored using a chemiluminescent analyzer for NOx, NDIR for C0/C02 
and a heated Flame Ionization Detector for HC as described by Nehmer [29]. 

The full-dilution tunnel was designed following EPA recommendations [29], 
and a schematic diagram of the apparatus is presented in Fig. 19. Sample 
emissions measurement results are presented in Fig. 20. In this case the 
engine was operated using the standard oil-test engine pump-line fuel 
injector for different runs and the results demonstrate the good repeatability 
of the emissions data. 



Fig. 18 Measured and KIVA predicted gas temperatures versus crank angle in 
the Caterpillar engine at a point 12 mm below the head, 25 mm from the axis. 
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Fig. 19 Schematic diagram of full-dilution tunnel for particulates 
measurements. 


CATERPILLAR 1Y540 

1900 rev/min, BSFC 0.348 lb/bhp-hr, BMEP 9.5 bar 



Run 

Fig. 20 Particulate and NOx measurements at 1900 rev/min showing 
repeatability of the emissions measurements for different runs. 


The engine experimental data are currently being used to test KTVA’s soot 
and emissions model predictions using the Nippondenso fuel injector with 
variable rate-shapes and injection timings, including split injections [29]. 
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EFFECT OF .INJECTION RATE AND SPLIT INJECTIONS - Electronically 
controlled diesel fuel injection equipment has been under active 
development for the past few years. Recent developments in microprocessor 
technology have created opportunities for better control of fuel delivery an 
timing than has hitherto been possible. This improved flexibility over the 
injection parameters offers the possibility of major improvements in diesel 
engine performance [30]. Improved injection timing flexibility and fast 
closing times are factors that are known to impact engine emissions [31, 32]. 
However, the improved flexibility also introduces more configuration 
possibilities, and this complicates the task of selecting optimum injection 
parameters. An objective of this research program is to apply advanced 
modeling techniques to provide guidance for the selection of fuel injection 
and combustion parameters. This will help in the evolution of more fully 
optimized, clean and efficient diesel engines. 

To help with the model validation effort, experiments are being conducted 
using the Nippondenso U2 injector system. A schematic diagram of the 
injector is presented in Fig. 21. The injector features a three-way solenoid 
valve with an interchangeable check valve for flow control. The initial 
injection rate can be adjusted by using orifices with different diameter chec 
valves. This is shown in Fig. 22 where check valve orifices of 0.6, 0.3 and 0.2 
mm diameter were used to produce injections with rise rates of 2, 7 and 16 
crank angle degrees to maximum open (at 1600 rev/min). The injector was 
operated at 90 MPa (13,230 psia) injection pressure and the nozzle tip features 
6 holes with hole diameter of 0.260 mm. The sprays are oriented with the 
spray axis at an angle of 27.5° to the head. The experimental conditions of the 
study are given in Table 8 [29]. 

Table 8. Experimental conditions for Ni ppondenso injector study 


Engine speed 
Overall equivalence ratio 
Air flow rate 
Fuel flow rate 
Intake tank pressure 
Intake tank temperature 
Exhaust back pressure 
Injection starts (nominal)* 
Injection ends (nominal) 

3 injection rates 

4 split injections 

2 different timings between 


1600 rev/min 

0.45 

63g/s 

2g/s 

1.82 atm 

36C 

1.57 atm 
10° BTDC 
15° ATDC 

( 2° , 7° , 16° to maximum open 
(0.6, 0.3, 0.2 mm check valve)) 
(10-90, 25-75, 50-50 and 75-25) 

(3 fi and 8° ) 


injections 


* Injection timing optimized for best BSFC at nominally constant fuel 
flow 
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Fig. 23 Measured effect of injection rate on engine cylinder pressure. 

The effect of variations in the injection rate on the engine performance is 
summarized in Fig. 23. As can be seen, the slowest injection rate case (0.2 mm 
check valve diameter) required the largest injection advance (injection was at 
20 degrees before TDC) and the BSFC was 194 g/kw-hr. The best fuel economy 
was obtained with the fastest rate, 0.6 mm check valve, where the injection 
was at 11 degrees BTDC and the BSFC was 183 g/kw-hr. 

The effect of injection rate on engine soot and NOx emissions is shown in 
Fig. 24. Comparing the data points labeled 'Base' (open circles) it is seen that 
the fast injection rate case, Base.6, with the 0.6 mm check valve, gave the 
lowest NOx and soot levels at about 0.1 and 2.4 g/hp-hr, respectively. These 
data are within the 1994 heavy duty truck standards, as is indicated in the 

figure. 

The soot (particulate) levels for the slower injection cases, Base.3 and Base.2, 
are similar to those obtained for the fast rate case, Base.6. However, the data 
indicates that the NOx levels increase with decreasing injection rate. The low 
NOx levels associated with the fast injection rate would be explained by the 
fact that this case also had the shortest ignition delay. 

The Soluble Organic Fraction (SOF) data are also given in Fig. 24. It can be 
seen that the SOF is within 10-15% of the particulate mass, indicating that the 
test engine is representative of a modem low SOF engine. 
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Fig. 24 Measured effect of injection rate (Base.2, .3 and .6) and split injections 
(10-90, 25-75, 50-50, 75-25% split) on Particulate and NOx emissions (see text). 

The effect of split injections on engine performance has also been studied. As 
indicated in Table 8, split injections were considered with the first injection 
containing 10, 25, 50 and 75% of the injected mass (and the second injection 
containing 90, 75, 50 and 25% of the injected mass, respectively). Time 
periods between the two injections of 3 and 8 crank angle degrees were 
considered. The corresponding injector needle lift and fuel injection rate data 
are presented in Fig. 25. As can be seen, the injector affords remarkable 
control over the details of the injection. 

The corresponding measured cylinder pressure data are shown in Fig. 26. 
Here it is seen that the injection details have a large influence on the 
combustion process. The injection timing was earlier, and the BSFC 
increased, as the amount of fuel in the first injection pulse was decreased for 
both of the spadngs between injections considered in the experiments. Figure 
24 shows that split injections also influence the engine particulate and NOx 
emissions significantly. In particular, the soot levels increased as the amount 
of fuel in the first injection decreased. However, the NOx emissions 
decreased as the amount of fuel in the first injection was decreased. The effect 
of the spacing between injections on the emissions is seen to be small. 
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Fig. 26 Measured cylinder pressure data for split injection cases. Left - 8 degree 
crank angle spacing between split injections. Right - 3 degree spacing. 


FLOW VELOCITY MEASUREMENT - Particle Image Velocimetry (PIV) is 
being used to make gas velocity and turbulence measurements in the 
(initially) motored engine. PIV provides velocity and turbulence data in a 2-D 
plane during one realization of the flow (not just at one point as with 
conventional laser Doppler velocimetry (LDV)). This later feature removes 
uncertainties associated with the cyclic variability of engine flows that are 
commonly encountered with single point measurements. The results will be 
used to validate the KIVA computations of the intake process. 

Similar engine modifications to those used for the schlieren temperature and 
diffusivity measurements will be used for the PIV study. The measurements 
of gas velocity are made in a plane within the piston bowl using the 
configuration shown in Fig. 27. A laser sheet is reflected off the central piston 
cusp which is coated with a reflective coating. A TSI particle seeder is used to 
supply TiC>2 particles to the intake flow to act as flow markers. A second Q- 
switch driver (Marx bank) has been installed on the Nd:YAG laser, and the 
laser is Q-switched twice during the fluorescence lifetime of the flash lamp in 
a preset time sequence. This produces two images of the Ti02 particles in the 
plane of the laser sheet on the 35 mm film. The flow velocities are then 
determined knowing the time-spacing of the images, and by measuring the 
distance moved by the particles using image analysis techniques. 

Preliminary PIV photographic data has already been obtained in a turbulent 
air jet, and data reduction software has been used to construct the details of 
the flow velocity field shown in Fig. 28. The details of the turbulent flow 
structures are clearly visible in the results. 
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Fig. 27 Schematic diagram of optical layout for in-cylinder PIV. 



Fig. 28 Velocity vectors in a 2-D plane in a turbulent air jet using PIV. 
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FUTURE WORK 


The research work for the next reporting period will address the entries in 
Table 1 (page 6) scheduled for year 4: intake flow, vaporization, and soot and 
radiation. In addition, work will continue on the combustion model and the 
engine experiments. A major focus will be the application of the model to 
the conditions of the engine experiments. Details of the individual tasks are 
described briefly below. 

Intake flow modeling 

Work on this project is continuing along several lines. First, more detailed 
experimental data is being obtained to help with validation. Experimental 
flow visualization, possibly with PIV data, will be used to compare large scale 
structures. This will help address issues pertaining to the turbulence 
modeling and grid resolution. Secondly, more unsteady calculations are 
being made and the flow structure is being examined in more detail. Particle 
traces and vertical structure of turbulence and other quantities are being 
obtained. Finally, since the unsteady intake calculations are expensive (over 
40 hours of Cray Y-MP time per run), efforts are under way to use the intake 
flow fields as initial conditions for other calculations. This will provide 
much better initial conditions for spray and combustion modeling studies. 

Multicomponent Fuel Vaporization 

The original KIVA vaporization model is a low pressure model that is based 
on the following assumptions : a. spherical symmetry, b. single component 
fuel, c. uniform drop temperature, d. constant fuel density and, e. unity Lewis 
number. Assumptions d. and e. have been removed in the present effort. 
The preparatory phases of modeling multi-component fuels (assumptions b. 
and c.) has been completed. The spherical symmetry assumption was 
considered to be reasonable for stable, non-breaking drops. 

Work is continuing on implementing a multi-component fuel vaporization 
model into KIVA that is based on the works of Refs. [33, 34]. 

Soot and Radiation 

The performance of the Surovkin and Hiroyasu soot models in engine 
computations will continue to be assessed. It is planned that the radiation 
model of Chang and Rhee [35] will be adapted for use in this study. This 
model assumes the Rayleigh-limit for the soot emissivity and a band model 
for the gas emissivity. More complicated radiation models that require 
solution of the radiation transport equation and that account for scattering by 
the spray drops are available [36]. However, in view of the uncertainties 
about soot concentrations, this additional complexity is not considered to be 
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warranted at_the present time. Much work is being done in this area and it 
will be monitored for its appropriateness for inclusion in KTVA. 

Combustion 

As discussed earlier, the present results indicate that the ignition process is 
well represented by the Shell multi-step kinetics model. However, the results 
indicate that the single step Arrhenius combustion model is inadequate for 
the main stage of combustion after ignition. To address this shortcoming the 
effect of turbulence will be included in the combustion model. This will be 
accomplished using a Magnussen-Hjertager type model [37] and also the 
characteristic time model of Kong et al. [12]. 

Experiments 

Engine experiments will be conducted to assess the effect of injection 
characteristics (injection rate and split injections) on engine performance 
over a wide range of engine operating conditions. It is planned to include 
combustion visualization experiments to allow the location and extent of 
combustion within the combustion chamber to be compared with the model 
predictions. These experiments will be conducted using a modified version 
of the engine with one of the exhaust valves replaced by a window for optical 
access. The PIV measurements of in-cylinder turbulence and velocity 
distributions will be conducted, and the results compared with the KIVA 
intake flow predictions. 

Model Application 

The model will be used to study the influence of engine operating conditions 
and fuel injection parameters (injection pressure, timing and rate shape - 
including pilot injection) on in-cylinder fuel-air distributions, combustion 
and emissions. The model predictions will be validated by comparison with 
the experimental cylinder pressure and other data. 

The rapid pressure rise that can result from premixed combustion of fuel that 
is vaporized during the ignition delay period has important implications on 
engine noise and structural loading. Over-mixed vaporized fuel is 
responsible for the increase in HC emissions found with an increased ignition 
delay. The influence of the fuel injection parameters on fuel vaporization 
and mixing during the ignition delay and the early stages of combustion will 
be studied using the model. Particulates are thought to be formed when 
under-mixed regions are exposed to high temperatures. The influence of the 
fuel injection parameters and in-cylinder flow details on mixing and 
combustion during the main burning phase will also be assessed. 
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SUMMARY AND CONCLUSIONS 

The three-dimensional KIVA code is being modified to include state-of-the- 
art submodels for diesel engine flow and combustion. Improved and/ or new 
submodels for wall heat transfer with unsteadiness and compressibility, 
laminar-turbulent characteristic time combustion with unburned HC and 
Zeldo'vich NOx, spray/ wall impingement with rebounding and sliding drops 
have already been implemented, as described in a previous annual report [5]. 

The present report describes progress on implementing improved models for 
spray drop drag and breakup, diesel ignition/combustion and soot formation. 
The results show that the effect of drop distortion on drop drag needs to be 
included for accurate spray computations. Also, the drop breakup results 
indicate that a surface wave breakup model predicts atomization drop sizes 
more accurately than the standard KIVA oscillating drop breakup (TAB) 
model. The combustion results indicate that the Shell multistep kinetics 
model gives accurate predictions of the diesel ignition process. However, 
further work is needed to assess the effects of turbulence on combustion 
during the later stages of combustion. The preliminary soot and NOx model 
results presented here are encouraging, but need further testing by 
comparison with engine experiments. 

A block structured version of KIVA together with an advanced grid 
generation scheme has been implemented. These codes have been applied to 
model a realistic engine geometry with steady and moving intake valves. 
Predicted steady flows agree reasonably with steady-state experimental data, 
and the moving valve results are encouraging. The accuracy of these results 
will be tested by comparison with measurements. 

An engine test facility has been constructed to allow validation of the in- 
cylinder flow and soot and NOx emission predictions. A single-cylinder, 
modem, well characterized research engine that meets the 1994 particulate 
and NOx emission standards has been setup. The engine's SOF levels were 
between 10-15 % of the soot mass which is also representative of the values 
found in modern low emission engines. 

The engine has been used to study the effect of injection rate and split 
injection on engine emissions. It was found that a fast rate of injection gives 
best overall engine performance with respect to load, soot and NOx 
emissions. In addition, the results show that a split injection has a large effect 
on the rate of cylinder pressure rise and the NOx emission is strongly effected 
by rate of cylinder pressure rise. Increased soot levels are obtained when a 
small first injection quantity (pilot injection) is used, possibly due to poor 
spray development, and a shorter interval between split injections improves 
BSFC. However, the effect of the interval between split injections on soot and 
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NOx is mixed. These results form a data base for future KIVA validation 
studies. 

An optoacoustic measurement technique has been developed for in-cylinder 
gas temperature and turbulent diffusivity measurements. Local gas 
temperature measurements in the motored diesel engine were in good 
agreement with KIVA model results. Other results in quiescent and 
turbulent jets indicate that the gas temperature measurements are accurate to 
within 3%. Work is in progress on improved optical arrangements for 
engine turbulent diffusivity measurements. 
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APPENDIX 1 KIVA USERS GROUP ACTIVITY 

The KIVA computer code and its predecessor codes were written at the Los 
Alamos National Laboratories, New Mexico [1,2]. These codes have seen 
increasing use by university researchers, government laboratories and engine 
industries in the U.S. and abroad since the early 1970s. A U.S. KIVA Users 
group was formed during the DOE Diesel Working Group Meeting at the 
University of Wisconsin-Madison in the Fall of 1989 . The purpose of this 
Users Group is to promote the use of KTVA and to facilitate exchange of 
information among KIVA users. Since that time KIVA Users Groups have 
also been formed in Europe and Japan. 

The U.S. KIVA Users Group currently numbers about 80 organizations. We 
have served as the editors of the KIVA Users Newsletter and seven issues 
have already been distributed to the membership (copies of newsletters 5, 6 
and 7 are attached in this Appendix). We also co-organized the First and 
Second International KIVA Users Group Meeting with CRAY Research. 
These meetings took place in February 1991 and 1992 in Detroit, MI. The 
meetings were attended by 75 to 80 registrants each. The Third International 
KIVA Users Group Meeting will take place on February 28, 1993 in Detroit. 
Details will be released in the 8th KIVA Users Newsletter. 
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0.4. For more information about this study contact — * L 12 — «b l 12 capabilities. For more information contact Reza 

the authors; Alex Uu and Rolf Reitz, Department of Fig. 3 Measured and predicted (parent) drop Taghavi, Industry, Science and Technology 
Mechanical Engineering University of Wisconsin, trajectories. The drops enter from the upper right Department, Cray Research Park, 655E Lone Oak 
Madison, Wl 53706. Phone: (608) 262 0145. and the measured data is indicated by the circles. Drive. Eagan, MN 55121. Phone: (612) 683-3643. 





APPENDIX 2 LIST OF RELATED PUBLICATIONS AND ABSTRACTS 


C. J. Rutland, C. M. Pieper and R. Hessel, "Intake and Cylinder How Modeling 
with a Dual- Valve Port," SAE Congress and Exposition, 1993 

Intake port and cylinder flow have been modeled for a dual intake valve diesel engine. A block 
structured grid was used to represent the complex geometry of the intake port, valves, and 
cylinder. The calculations were made using a pre-release version of the KIVA-3 code developed 
at Los Alamos National Laboratories. Both steady flow-bench and unsteady intake 
calculations were made. 

In the flow bench configuration, the valves were stationary in a fully open position and pressure 
boundary conditions were implemented at the domain inlet and outlet. Detailed structure of the 
in-cylinder flow field set up by the intake flow was studied. Three dimensional particle trace 
streamlines reveal a complex flow structure that is not readily described by global parameters 
such as swirl or tumble. Streamlines constrained to lie in planes normal to the cylinder axis 
show dual vortical structures, which originated at the valves, merging into a single structure 
downstream. Initial comparisons of the computational results to data from steady flow-bench 
experiments are encouraging. 

In the unsteady intake calculations, dual moving valves and a piston with a bowl were 
implemented. Calculations start prior to intake valve opening and end just prior to fuel 
injection. Velocity vector plots show the characteristic jet like flows from the valve curtain and 
vortical flow structures in the combustion chamber. Volume averaged turbulent kinetic energy 
and length, and swirl ratio magnitude increase during intake and decay during compression. The 
swirl structure is examined and compared to mean values by averaging over horizontal planes. 


R.P. Giangregorio, Y. Zhu, R.D. Reitz, "Application of Schlieren Optical 
Techniques for the Measurement of Gas Temperature and Turbulent 
Diffusivity in a Diesel Engine," SAE Congress and Exposition, 1993 


A new technique which is based on optoacoustic phenomena has been developed for measuring 
in-cylinder gas temperature and turbulent diffusivity. In the experiments, a high energy 
Nd:YAG pulsed laser beam was focused to cause local ionization of air at a point in the 
combustion chamber. This initiates a shock wave and creates a hot spot. The local temperature 
and turbulent diffusivity are determined by monitoring the shock propagation and the hot spot 
growth, respectively, with a schlieren photography system. 

In order to assess the validity and accuracy of the measurements, the technique was also 
applied to a turbulent jet. The temperature measurements were found to be accurate to within 
3%. Results from the turbulent jet measurements also showed that the growth rate of the hot 
spot diameter can be used to estimate the turbulent diffusivity. 

In-cylinder gas temperature measurements were made in a motored single cylinder Caterpillar 
diesel engine, modified for optical access. The measured results were in good agreement with 
those from multidimensional calculations using the KIVA model. However, the presence of 
density gradients in the gas due to the compression process was found to deteriorate the 
schlieren images of the hot spot in the engine, and alternative optical arrangements are being 
explored for engine turbulent diffusivity measurements. 
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A.B. Liu, D. Mather, and R. D. Reitz, "Modeling the Effects of Drop Drag and 
Breakup on Fuel Sprays," SAE Congress and Exposition, 1993 

Spray models have been evaluated using experimentally measured trajectories and drop sizes of 
single drops injected into a high relative velocity gas flow. The computations were made using 
a modified version of the KIVA-2 code. It was found that the drop drag coefficient and the 
drop breakup time model constant had to be adjusted in order to match the measurements. 
Based on these findings, a new drop drag submodel is proposed in which the drop drag 
coefficient changes dynamically with the flow conditions. The model accounts for the effects of 
drop distortion and oscillation due to the relative motion between the drop and the gas. The 
value of the drag coefficient varies between the two limits of that of a rigid sphere (no 
distortion) and that of a disk (maximum distortion). The modified model was also applied to 
diesel sprays. The results show that the spray tip penetration is relatively insensitive to the 
value used for the drop drag coefficient. However, the distribution of drop sizes within sprays 
is influenced by drop drag. This is due to the fact that changes in drop drag produce changes in 
the drop-gas relative velocity. This, in turn, causes changes in the spray drop size through the 
drop breakup and coalescence processes. The changes occur in such a way that the net effect on 
the spray penetration is small over the tested ranges of conditions. These results emphasize 
that measurements of spray penetration are not sufficient to test and produce improved spray 
models. Instead, local measurements of drop size and velocity are needed to develop accurate 
spray models. 

R.D. Reitz, N. Ayoub, M. Gonzalez, R. Hessel, S. Kong, J. Lian, C. Pieper and 
C.J. Rutland, "Improvements in 3-D Modeling of Diesel Engine Intake Flow 
and Combustion," SAE Paper 921627, 1992 

A three-dimensional computer code (KIVA) is being modified to include state-of-the-art 
submodels for diesel engine flow and combustion: spray atomization, drop breakup/coalescence, 
multi-component fuel vaporization, spray/ wall interaction, ignition and combustion, wall heat 
transfer, unbumed HC and NOx formation, soot and radiation and the intake flow process. 
Improved and/or new submodels which have been completed are: wall heat transfer with 
unsteadiness and compressibility, laminar-turbulent characteristic time combustion with 
unbumed HC and Zeldo'vich NOx, and spray/wall impingement with rebounding and sliding 
drops. Results to date show that: adding the effects of unsteadiness and compressibility 
improves the accuracy of heat transfer predictions; spray drop rebound can occur from walls at 
low impingement velocities (e.g., in cold-starting); larger spray drops are formed at the nozzle 
due to the influence of vaporization on the atomization process; a laminar-and-turbulent 
characteristic time combustion model has the flexibility to match measured engine combustion 
data over a wide range of operating conditions; and, finally, the characteristic time combustion 
model can also be extended to allow predictions of ignition. The accuracy of the predictions is 
being assessed by comparisons with available measurements. Additional supporting 
experiments are also described briefly. To date, comparisons have been made with measured 
engine cylinder pressure and heat flux data for homogeneous charge, spark-ignited and 
compression-ignited engines, and also initial comparisons for diesel engines. The model results 
are in good agreement with the experiments. 
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S.-C. Kong, and R. D. Reitz, "Multidimensional Modeling of Diesel Ignition 
and Combustion Using A Multistep Kinetics Model,” ASME Internal 
Combustion Engine Symposium, Energy-sources Technology Conference and 
Exhibition, January 31-February 4, 1993, Houston, TX. 

Ignition and combustion mechanisms in diesel engines were studied using the KIVA code, with 
modifications to the combustion, heat transfer, crevice flow and spray models. A lammar-and- 
turbulent characteristic-time combustion model that has been used successfully for spark- 
ignited engine studies, was extended to allow predictions of ignition and combustion in diesel 
engines. A more accurate prediction of ignition delay was achieved by using a multistep 
chemical kinetics model. The Shell model was implemented for this purpose and was found to 
be capable of predicting successfully the autoignition of homogeneous mixtures in a rapid 
compression machine and diesel spray ignition under engine conditions. The physical 
significance of the model parameters Is discussed and the sensitivity of results to the model 
constants is assessed. The ignition kinetics model was also applied to simulate the ignition 
process in a Cummins diesel engine. The post-ignition combustion was simulated using both a 
single-step Arrhenius kinetics model and also the characteristic-time model to account for the 
energy release during the mixing-controlled combustion phase. The present model differs from 
that used in earlier multidimensional computations of diesel ignition in that it also includes 
state-of-the-art turbulence and spray atomization models. In addition, in this study the model 
predictions are compared to engine data. It is found that good levels of agreement with the 
experimental data are obtained using the multistep chemical kinetics model for diesel ignition 
modeling. However, further study is needed of the effects of turbulent mixing on post-ignition 
combustion. 


56 


APPENDIX 3 FORTRAN SUBROUTINES 

The FORTRAN subroutine given in this appendix has been tested for diesel 
sprays with combustion. The subroutine needs to be compiled and linked 
together with the standard KTVA subroutines. 

Subroutine chem (replaces KIVA subroutine chem) 

This subroutine is based on the Ref. [17] Theobald, M.A., and Cheng, W.K., "A 
Numerical Study of Diesel Ignition," ASME Paper 87-FE-2, 1987 and should be 
used with Subroutine Chempm [see Ref. 5]. 

subroutine chem 
c based on 

Theobald, M.A., and Cheng, W.K., "A Numerical Study of Diesel 
Ignition," ASME Paper 87-FE-2, 1987. 

csong modified for the use of SHELL autoignition model combined with 
c charateristic-time combustion model which is ' chemprn . f ' 

c apply in combustion bomb spray ignition and Cummins NH engine 

c 

c (1) species order: 1: fuel 2: 02 3: N2 4: C02 5: H20 

c 6: R 7: B 8: Q 9: CO 10: H2 

c (2) in 'itape' set the enthalpy of formation of new species, i.e., 
c R, B, Q, equal to zero. 

c (3) in 'rinput.f' change the enthalpy of new species, i.e., 
c R, B, Q, equal to zero. 

c (4) switch the low-high-temperature chemistry at TCHOP=1000 K 
c (5) the high temperature chemistry refers to the use of the 

c characteristic-time combustion model. 

c (6) species order in previous characteristic-time combustion model 
c has to be changed accordingly. 

c 

c 

c chem is called by: kiva 

c 

c chem calls the following subroutines and functions: chemprn 

c 

c 

c 

include 1 comd . com ' 

CC&MDN/CHUBS/ICSLOG, ICSLC 
common/hrr/tothrr 
dimension domega (lnrk) 
real kf , kb 

DATA X1,X3,X4 /l. , 0., -1 . /, Yl,Y3,Y4/0. , 0. , 0.35/ 
cl2h26 dodecane 2d spray 

C DATA FUELN, FUELM, FUELQ, COF, C02F , PSTOIC, CRATIO, QFUEL, RMN2 
c &/12., 13., 1.92308, . 61846, . 30462, 1 . 11385, . 67, 3. 9363E12, 3 . 60714/ 
cl4h30 tetradecane 3d Cummins 

DATA FUELN, FUELM, FUELQ, COF, C02F, PSTOIC, CRATIO, QFUEL, RMN2 
&/14. ,15. ,1.93333, .62533, .30800,1.12100, .67,3. 9363E12, 4. 10714/ 

DATA AF01,EF1,AF02,EF2,AF03,EF3,AF04,EF4 
&/7 . 3E-4, -7 . 55E3, 180 . , -3 . 53E3, 1 . 47, 5. 04E3, 1 . 30e6, 1 . 51E4/ 

DATA TFREEZ , RTFREZ , TINHIB, CQINH 
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& /950. , 1 . 052631E-3, 1000 . , 1 . E-7/ 

DATA AP1,EP1, AP2,EP2,AP3,EP3 
& /1.E12 , 0., 1 .Ell, 7.56E3, 1.E13, 428./ 

DATA IGLOG/0/, ICSLOG/O/, IGLOGO/O/, IGIXX32/0/ 
DATA TCHOP /1000./ 

DATA XSCHEM /1.00/, ICSLIM/200/ 



c 

tchem=l . 0e-10 

icslog=0 

tothr=0 . 

do 100 k=l,nz 

i4b= (k-1) *nxpnyp 

do 90 j=l, ny 

i4=i4b+ ( j-1) *nxp+l 

do 80 i=l,nx 

if (f (14) . eq.0. ) go to 80 

if (spd(i4, 1) .le.0.) spd(i4,l)=l.e-20 

if (spd(i4, 2) . le.0.) spd(i4,2)=l.e-20 

IF (TEMP (14) .LT.TCUT)GO TO 80 

SISTOC—SIE (14) 

C DO NOT NEED TO SUBCYCLE IF [B] TOO LOW FOR BRANCHING 
IF (SPD (14, 7) *RMW(7) .LE. 1 .E-15) GO TO 9000 

C DO NOT SUBCYCLE IF CELL TOO HOT FOR AUTOIGNITION EQUATIONS 
IF (TEMP (14) .GT. TCHOP) GO TO 9000 
RTI JK=1 . /TEMP (14) 

IF (TEMP (14) .GT.TFREEZ) RTIJK=RTFREZ 
TKBEST=CF (4) *EXP (-EF (4) *RTIJK) 

C TOTAL NET B RATE CONSIDERED TO SET SUBCYCLING 
CONB= SPD (14,7) *RMW (7) 

CONR= SPD (14, 6) *RMW(6) 

CONQQ= SPD (14,8) *RMW (8) 

CONRH= SPD (14,1) *RMW (1) 

CONOX= SPD (14, 2) *RMW (2) 

P1K= API* EXP (-EP1*RTIJK) 

P2K= AP2*EXP (~EP2*RTI JK) 

P3K= AP3*EXP(-EP3*RTIJK) 

PTK=1 . / (1 . / (CONOX*PlK) +1 . /P2K + 1 . / (CONRH*P3K) ) 
F01=AF01*EXP (-EF1*RTI JK) 

F02=AF02*EXP (-EF2*RTIJK) 

F1L=F01* (CONOX**Xl) * (CONRH**Yl) 


TKBEST=ABS ( (-TKBEST*CONB+F1L*PTK*CONR+F2P*CONR*CONQQ) /CONB) 

ICS=IFIX (DT*TKBEST/XSCHEM) 

if (ics.le.l) ics=l 

IF (ICS . GT . ICSLIM) ICS=ICSLIM 

DTC=DT/ (FLOAT (ICS) ) 


GO TO 9010 
9000 ICS=1 
DTC=PT 

9010 CONTINUE 

IF (ICS . GT . ICSLOG) ICSLC=I4 
IF (ICS . GT . ICSLOG) ICSLOG=ICS 

+++ SET EXPONENTS AND STOICH. COEF. FOR SHELL .... 

DO NOT MOVE THESE TO A DATA STATEMENT AS THEY ARE 
OVERWRITTEN IN RINPUT 
FBM (4,2) =C02F 
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FBM(9,2)=COF 
FBMAM(4,2)=C02F 
FBMAM (9,2) =COF 
FBM (3, 5) =RMN2 
FBM(3,6)=2.*RMN2 
FBMAM (3, 5)=RMN2 
FBMAM (3, 6) =2 . *RMN2 

c +++ 

DO 75 IC=1,ICS 

C SECOND STAGE IGNITION DETECTOR 

IF ((TEMP (14) .GT. 1000 . ) .AND. (IGLOGO.EQ.O) )GO TO 7002 
GO TO 7003 

7002 WRITE ( *,7012) T+ (IC-1) *DTC, 14 
IGLOG0=l 

7003 IF ( (TEMP (14) .GT.1100. ) .AND. (IGLOG.EQ.O) )GO TO 7000 
GO TO 7004 

7000 WRITE ( *,7010) T+ (IC-1) *DTC, 14 

7010 FORMAT (5X, *1100 DEGREE K IGNITION AT T=',E10.4,' AT CELL ’,15) 
IGLOG=l 

7004 IF ( (TEMP (14) . GT. 2000 . ) . AND . (IGLOG2 .EQ. 0) ) GO TO 7006 
GO TO 7020 

7006 WRITE ( *,7014) T+ (IC-1) *DTC, 14 
IGLOG2=l 

7020 CONTINUE 

7012 FORMAT (5X, '1000 DEGREE K IGNITION AT T=',E10.4,' AT CELL ',15) 

7014 FORMAT (5X, '2000 DEGREE K IGNITION AT T=',E10.4,' AT CELL ’,15) 
ti jk=temp(i4) 
if (tijk.lt. tcut) go to 80 
rti jk=l. /ti jk 

FREEZE REACTION RATES OF KNOCK MODEL TO 950 K IF TEMP HIGHER. 

(NOTE PRE -EXPONENTIAL TEMP FACTOR NOT FROZEN, BUT NEVER USED) 
IF (TI JK. GT . TFREEZ) RTI JK=RTFREZ 
CONQQ=SPD (14, 8) *RMW (8) 

C SHUTDOWN AUTOIGNITION EQUATIONS IF CELL ALREADY HOT. 

IF (TIJK.GT.TCHOP) GO TO 900 
CONOX=SPD (14, 2) *RMW (2) 

CONRH=SPD (14, 1) *RMW (1) 

C JUMP OUT OF LOOP IF NO FUEL OR 02 AVAILABLE 

IF (CONRH.LT. l.E-10) GO TO 76 
IF (CONOX.LT. l.E-10) GO TO 76 
F01=AF01*EXP (-EF1*RTIJK) 

F02=AF02*EXP (-EF2*RTIJK) 

F03=AF03*EXP (-EF3*RTIJK) 

F04=AF04*EXP (-EF4*RTIJK) 

P1K= AP1*EXP (-EP1*RTI JK) 

P2K= AP2*EXP(-EP2*RTIJK) 

P3K= AP3*EXP(-EP3*RTIJK) 

F1L=F01* (CONOX**Xl) * (CONRH**Yl) 

F3L=F03* (CONOX**X3) * (CONRH**Y3) 

F4L=F04* (CONOX**X4) * (CONRH**Y4) 

PTK=1 . / (1 . / (CONOX*PlK) +1 . /P2K+1 . / (CONRH*P3K) ) 

C ABOVE FORM FOR PTK IS CORRECTED VERSION 11/27/85 
F1P=F1L*PTK 
F2P=F02*PTK 
F3P=F3L*PTK 
F4P=F4L*PTK 

C SPECIAL TEMP/CONC. DEPENDENT STOICHIOMETRY FACTORS 
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EQLL= (F1L*MW (7) +F4L*MW(8) ) / (MW(1) /FUELM+PST0IC*MW(2) ) 
ASP1= (EQLL+1 . ) /FUELM 
ASP2= (EQLL+1 . ) +PSTOIC 

C SET STOICH. FACTORS IN PROPAGATION EQUATION NRK=2 
FAM(1, 2) = ASP1 
FAM(2, 2) = ASP2 
FBM (7,2) = F1L 
FBM (8,2) = F4L 
FBMAM ( 1 , 2) = -ASP1 
FBMAM (2,2) = -ASP 2 
FBMAM (7, 2) = F1L 
FBMAM(8, 2) = F4L 
DO 70 IR=1 , NRK-1 
rp=l . 
pp=l. 

ne=nelem(ir) 

do 20 kk=l,ne 

isp=cm (kk, ir) 

rom=spd(i4, isp) *rmw(isp) 

if (am(isp,ir) .eq.O) go to 10 

if (rom.le.O. ) rp=0. 

if (rom.gt.O.) rp=rp*rom**ae (isp, ir) 

10 if (bm(isp,ir) .eq.O) go to 20 
if (rom.le. 0. ) pp=0. 
if (rom.gt.O.) pp=pp*rom**be (isp, ir) 

20 continue 
kb=0 . 
kf=0. 
teback=l . 
teford=l. 
ekback=l . 
ekf ord=l . 

30 if (cf (ir) . le . 0 . ) go to 40 
c +++ 

c +++ forward reaction coefficient 
c +++ 

if (ef (ir) .ne.0.) ekford=exp (-ef (ir) *rti jk) 
if (zetaf (ir) .ne.0. ) teford=ti jk**zetaf (ir) 
kf=cf (ir) *teford*ekford 
IF ( IR . EQ . 2 ) KF=PTK 
IF (IR.EQ. 3) KF=F2P 
IF ( IR . EQ . 5 ) KF=F3P 

c +++ 

c +++ if any rate coefficients cannot be put in standard 
c +++ form, code them by hand and put them here 
c +++ 

c +++ find the reference species (the one in greatest danger 
c +++ of being driven negative) 
c +++ 

40 omeg=kf*rp-kb*pp 
rmin=0. 

if (omeg.eq.0. ) go to 70 
do 50 kk=l , ne 
isp=cm(kk, ir) 

if (spd(i4, isp) . le . 0 . ) go to 50 
rom=omeg*fbmam(isp, ir) *mw(isp) /spd(i4, isp) 
if (rom.ge. 0. ) go to 50 
if (rom.lt .rmin) iref=isp 
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rrain=min (rmin, rom) 

50 continue 

rom=spd(i4, iref ) *rmw (iref ) 
flam=fam(iref , ir ) 
f lbm=fbm (iref, ir) 
ctop=f lam*kb*pp + flbm*kf*rp 
cbot=f lam*kf *rp + flbm*kb*pp 

domega (ir) =rom*dtc* (ctop-cbot) / ( (rom+dtc*cbot) * (f lbm-f lam) ) 
do 60 isp=l,nsp 

spd(i4, isp) =spd(i4, isp) +mw (isp) *fbmam(isp, ir) *domega (ir) 

60 continue 

HEAT RELEASE ONLY ON PROPAGATION STEP SCALED BY 
REACTION CONSUMPTION. 

DECHEM = 0. 

IF (IR.EQ.2) DECHEM=QFUEL* (EQLL+1 . ) *DOMEGA(2) /RO(I4) 
sie (i4) =sie ( i 4 ) +dechem 
70 continue 
900 CONTINUE 

C DISALLOW HEAT RELEASE IF T.LT.TINHIB [NEW] 

IF (TI JK. LT. TINHIB) GO TO 970 
csong switch to char, time model 071792 
call chemprn(i4) 
spd(i4, 6)=l.e-20 
spd(i4, 7)=l.e-20 
spd (i 4 , 8 ) =1 . e-20 
go to 970 

C MAIN HEAT RELEASE IR=NRK. LAST KINETIC EQN 
IR=NRK 

C RESET TEMP TO UNFREEZE IF PREVIOUSLY FROZEN 

RTI JK=1 . /TI JK 
RP=1. 

PP=1. 

NE=NELEM(IR) 

DO 920 KK=1, NE 
ISP=CM(KK, IR) 

ROM=SPD (14, isp) *RMW (ISP) 

IF (AM(ISP,IR) .EQ.0) GO TO 910 
IF (ROM.LE . 0 . ) RP=0 . 

IF (ROM. GT. 0 . ) RP=RP*ROM**AE (ISP, IR) 

910 IF (BM(ISP,IR) .EQ.0) GOTO 920 
IF (ROM.LE. 0 . ) PP=0. 

IF (ROM. GT . 0 . ) PP=PP*ROM**BE (ISP, IR) 

920 CONTINUE 
KB=0 . 

KF=0 . 

TEB=1 . 

TEF=1 . 

EKB=1. 

EKF=1 . 

IF (CB(IR) .LE.0. ) GO TO 930 
C BACKWARD REACTION COEF 

IF (EB (IR) .NE.0.) EKB=EXP (-EB (IR) *RTIJK) 

IF (ZETAB (IR) . NE . 0 . ) TEB=T I JK* * Z ETAB (IR) 

KB=CB (IR) *TEB*EKB 
930 IF (CF (IR) .LE. 0 . ) GO TO 940 
C FORWARD REACTION COEF 

IF (EF (IR) .NE.0. ) EKF=EXP (-EF (IR) *RTIJK) 
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IF ( ZETAF ( IR) . NE . 0 . ) TEF=TI JK* * ZETAF ( IR) 

KF=CF (IR) *TEF*EKF 
NONSTANDARD RATE COEF GO HERE 
FIND REF SPECIES 
940 OMEG=KF*RP-KB*PP 
RMIN=0. 

IF (OMEG .EQ. 0.) GO TO 970 
DO 950 KK=1, NE 
ISP=CM (KK, IR) 

IF (SPD (14, ISP) .LE . 0 . ) GO TO 950 
ROM=OMEG*FBMAM(ISP, IR) *MW (ISP) /SPD (14, ISP) 

IF (ROM.GE . 0 . ) GO TO 950 
IF (ROM.LT.RMIN) IREF=ISP 
RMIN=AMIN1 (RMIN, ROM) 

950 CONTINUE 

ROM=SPD (14, IREF) *BMW (IREF) 

FLAM=FAM (IREF , IR) 

FLBM=FBM ( IREF , IR) 

CTOP=FLAM*KB*PP+FLBM*KF*RP 

CBOT=FLAM*KF*RP+FLBM*KB*PP 

DOMEGA(IR) =ROM*DTC* (CTOP-CBOT) / ( (RQM+DTC*CBOT) * (FLBM-FLAM) ) 
DO 960 ISP=1,NSP 

SPD (14, ISP) =SPD (14 , ISP) +MW (ISP ) *FBMAM (ISP , IR) *DOMEGA (IR) 

960 CONTINUE 

DECHEM=QR(NRK)*DQMEGA(IR) /RO(I4) 


8000 


c +++ 

DECHK=ABS (DECHEM/SIE (14) ) 

SIE (14) =SIE (14) +DECHEM 
TCHEM=AMAXl (TCHEM, DECHK) 

970 CONTINUE 

IF (TEMP (14) .GT. 5000.) GOTO 
IF(ICS.EQ.l)GO TO 75 

8000 WRITE ( *,8110) T,NCYC,I, J, K, 14, TEMP (14) , SIE (14) , IT, CRANK 

8110 FORMAT ( * TEMPERATURE OVERFLOW AT T=',E12.5, 1 CYC' , 15, 1 1=' , 13, 
& *J=\I3, 'K=',I3, *14=', 15/' TEMP=',E12.5, ' SIE=',E12.5, 

& 'IT-' ,15, ' CRANK= * , 0PF7 . 2 ) 

75 CONTINUE 

C END OF SUBCYCLE LOOP ON CHEMISTRY 

76 CONTINUE 

DESIE=ABS ( (SIE (14) -SISTOC) /SIE (14) ) 

TCHEM=AMAX1 (TCHEM, DESIE) 

C END HEAT RELEASE COMPUTATION 

tothr=tothr+ (sie (14 ) -sistoc) *ro (i4) *vol (14 ) 
csong end of SHELL subcycling for current i4 at current time step. 

80 i4=i4+l 
90 continue 
100 continue 

tothrr=tothr*l . e-10/dt 

return 

end 
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